Finite-size version of the excitonic instability in graphene quantum dots 
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By a combination of Hartree-Fock simulations, exact diagonalization, and perturbative calcula- 
tions, we investigate the ground-state properties of disorder-free circular quantum dots formed in a 
graphene monolayer. Taking the reference chemical potential at the Dirac point, we study N < 15 
interacting particles, where the fine structure constant a parametrizes the Coulomb interaction. We 
explore three different models: (i) Sucher's positive projection ("no-pair") approach, (ii) a more gen- 
eral Hamiltonian conserving both N and the number of additional electron-hole pairs, and (iii) the 
full quantum electrodynamics (QED) problem, where only N is conserved. We find that electron- 
hole pair production is important for a > 1. This corresponds to a reconstruction of the filled Dirac 
sea and is a finite-size version of the bulk excitonic instability. We also address the effects of an 
orbital magnetic field. 

PACS numbers: 03.65.Pm, 73.22.Pr, 71.15.Bf, 73.21.La 
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I. INTRODUCTION 



Coulomb interaction effects in monolayer graphene 1 * 2 - 
are currently attracting a lot of attention; for a recent re- 
view, sec Ref.S From a theory point of view, this interest 
mainly stems from the possibility of realizing a strong- 
coupling version of QED in a readily accessible two- 
dimensional (2D) system. In fact, the (bare) fine struc- 



ture constant is rather large, a 



'/(Hkvf) ~ 2.2/k, 



with the effective substrate dielectric constant k and the 
Fermi velocity vp ~ c/300 sa 10 6 m/s. Retardation 
effects are irrelevant here, i.e., we effectively have 2D 
massless Dirac fcrmions interacting via the Coulomb po- 
tential. Similar physics can be expected for the surface 
state in 3D topological insulators,— but interactions are 
expected to be much weaker due to the large k in the rele- 
vant materials. In graphene, the situation away from the 
Dirac point (defined as zero of energy) can be reasonably 
well understood in terms of Fermi liquid theory,— £ but 
the picture is more complicated near the Dirac point. At 
a critical interaction strength a c , a semimetal- insulator 
transition is theoretically expected^ due to electron-hole 
proliferation. For a > a c , a finite gap corresponding 
to an excitonic insulator is formed and the ground state 
undergoes reconstruction. On the other hand, quantum 
critical behavior is expected 7 ^ as a precursor of the insta- 
bility for a < a c . Recent lattice quantum Monte Carlo 
simulations* found the critical value a c w 1.1 for an in- 
finitely extended ("bulk") graphene monolayer. Similar 
values were also obtained analytically from the dynamical 
polarization function approach 1 ^ and under the ladder 
approximation to the Bcthc-Salpctcr equation^ How- 
ever, so far no experimental signature of this excitonic 
instability has been reported. It has also been recog- 
nized that the excitonic instability for the bulk many- 
body problem is related to the simpler "supercritical" 
instability of the hydrogen problem in graphene j» where 
a corresponds to the (attractive) potential strength of 
the nucleus. Above a critical value for a, the nucleus cap- 
tures an electron to screen its positive charge below criti- 



cality, while at the same time a hole escapes to infinity in 
order to maintain charge neutrality. In atomic physics, 
essentially the same phenomenon should also take place 
for superheavy atoms.— The creation of an electron-hole 
pair thus also accompanies the supercritical instability. 
In the presence of an homogeneous orbital magnetic field 
B, the hole escape process is disturbed by the forma- 
tion of closed Landau orbits. For the bulk many-body 
problem, the resulting magnetic catalysis phenomenon 1 ^ 
implies a lowering of a c with increasing B. 

In this work, we study a finite-size version of the 
excitonic instability presumably realized in available 
graphene quantum dots. Quantum dots in conven- 
tional 2D systems have been studied extensively ) 14 ' 15 
and experimental results for lithographically prepared 
graphene dots were reported recentlyiiSr— Within the 
single-particle picture, theoretical proposals on how to 
model such a dot have been reviewed in Rcf. i We 
here adopt the probably simplest route by imposing the 
so-called "infinite-mass boundary condition," 21 ' 22 where 
no current is allowed to flow through the circle r = R 
defining the dot's boundary. While disorder limits the 
quality of the boundary in existing dots^ such a bound- 
ary condition captures at least their qualitative physics. 
Moreover, future experimental progress is likely to yield 
well-defined boundaries. 

We investigate the ground state of N interacting elec- 
trons in a closed circular graphene dot, where N parti- 
cles are added on top of the filled Dirac sea, i.e., rel- 
ative to the chemical potential /j, = 0. This prob- 
lem has been studied before within the Hartree-Fock 
(HF) approach i 2 ^— However, when going beyond ef- 
fective single-particle theory, one has to deal with the 
"Brown-Ravenhall disease," 27 ' 28 i.e., the possibility to 
excite electron-hole pairs with small energy by combin- 
ing a hole and an electron both very far away from the 
Fermi surface. While such processes are physically sup- 
pressed by the finite bandwidth, the infinitely deep filled 
Dirac sea present in the Dirac theory renders naive ap- 
proaches mathematically ill-defined. For a <C 1 and 
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when a gap separates electron and hole states, Sucher— 
showed that one can circumvent the Brown- Ravenhall 
problem by a suitable projection A + of the basic QED 
Hamiltonian H to a well-defined no-pair Hamiltonian 
H+ = A + HA + , where the filled Dirac sea is effectively 
treated as completely inert. The projection operator A+ 
eliminates negative-energy (hole) states from the single- 
particle Hilbert space. We study the validity of the no- 
pair approach for graphene dots and find that for a < 1, 
it is indeed meaningful, see also Ref. I29I On a quanti- 
tative level, however, it is accurate only for a <C 1. As 
also discussed by Sucher^ if one wishes to go beyond 
the positive projection scheme, a QED approach is in- 
dicated. The QED Hamiltonian H, see Eq. (fTTj) below, 
does not conserve the number N e h of electron-hole pairs. 
In fact, only the particle number N - defined as the im- 
balance of electron and hole numbers - is conserved, and 
a superposition of states with different N e h determines 
the ground state for strong interactions. Once electron- 
hole pairs proliferate, a reconstruction of the ground state 
takes place. We encounter this phenomenon for a > 1 in 
graphene dots, similar to the reported critical valued for 
the bulk excitonic instability. However, in our finite-size 
system this is a smooth crossover and not a phase tran- 
sition. We stress that our 'W-particle problem" defines 
N as the difference of electron and hole numbers, which 
allows for the excitation of an arbitrary number N e h of 
electron-hole pairs. For a — > 0, this definition reduces to 
having N electrons on top of the filled Dirac sea. 

The structure of the remainder of this paper is as fol- 
lows. In Sec. |H] we introduce the model and discuss 
the various theoretical approaches employed to find the 
ground state. An intermediate approach is to general- 
ize the no-pair approach (where N e h = 0) to allow for a 
fixed but finite number N e h of electron-hole pairs. The 
Hamiltonian is obtained from H by neglecting all 
terms that do not conserve N e h- A sufficient (but not 
necessary) condition for the breakdown of the no-pair 
Hamiltonian H+ arises when the ground-state energy of 
ilfix is lowered for some N e h > 0. We assess the validity 
of the no-pair scheme in Sec. IHII bv comparing to results 
obtained under iJfi x and from the QED Hamiltonian H. 
We perform these calculations using exact diagonaliza- 
tion (ED) for N = 2 and N = 3 particles in the dot. In 
Sec. lIVl we use H+ to carry out detailed HF calculations 
for up to N = 15 particles and relatively weak interac- 
tions, a < 1. We present results for the ground-state 
spin, valley polarization, and addition energy as function 
of N. Finally, in Sec. |V] we provide a discussion of our 
main results. 



II. MODEL AND THEORETICAL 
APPROACHES 

In this section, we describe the model employed in our 
study of the electronic properties of interacting graphene 
quantum dots. We will then turn to different theoretical 



approaches to obtain the ground-state properties. 

A. Single-particle problem 

It is well established that on low energy scales, 
quasiparticles in graphene are described by the Dirac 
Hamiltonian^ 

H = v F cr ■ (p + - A) + M(r)a z T z - ^ B s • B, (1) 

where p = —ih(d x , d y ) T and — e is the electronic charge. 
The Pauli matrices <x = (a x ,a y ) and a z refer to 
graphene's sublattice structure, while the Pauli matrix 
t z corresponds to the valley degree of freedom, i.e., to 
the two K points. A static vector potential A(r) [with 
r = (x,y) T ] allows for the inclusion of a constant orbital 
magnetic field B Zl where we choose the symmetric gauge, 
A = 7}B z (—y,x) T . Since we neglect spin-orbit couplings 
in Eq. ((TJ), spin Pauli matrices s = (s x ,s y ,s z ) only ap- 
pear in the Zeeman term. With fis denoting Bohr's mag- 
neton and putting the Lande factor to g e = 2^S this 
term couples to the full (homogeneous) magnetic field, 
B = (B x , By, B z ) with B = |B|. Switching to polar 
coordinates (r,(p), we consider a clean circular quantum 
dot in a graphene monolayer modelled by the well-known 
infinite- mass boundary condition^ where the mass M(r) 
in Eq. ([1} is zero for r < R but tends to +00 for r > R. 
This choice ensures that no current flows through the 
boundary at r = R. Eigcnstates can be classified by the 
conserved total angular momentum, J z = —ihd < p + Ha z /2, 
with eigenvalue hj for half-integer j — m + 1/2, meZ, 

While the eigenfunctions can be found in analyti- 
cal form even in the presence of the magnetic field^ 
the Coulomb interaction matrix elements are readily 
available^ only in the B = basis. We therefore first 
describe the solution for B = and later include the 
homogeneous magnetic field. Note that different valleys 
(t = ±) are decoupled and spin (s = ±) then simply 
yields a twofold degeneracy. For given (m, r, s), we first 
discuss the E > solutions to H <&( + ) = E& + \ where 
the spinor has the sublattice structure 

* (+) ^« = e ^(ie^L P (r))- ^ 
The infinite-mass boundary condition implie d 21 ' 22 

i>X, m {R) = Tll>2, m {R). (3) 

With the Bessel functions J m (kr) of the first kind, k = 
E/Hvp, and normalization constant A, the Dirac equa- 
tion for r < R is solved by the Ansatz 

ipi, m (r) = AJ m (kr), 4>2, m = AJ m+1 (kr). 

The quantization condition §3§ then determines the 
eigenenergies E a > with a = (n, m, t, s), 

J m {E a /A ) - TJ m+1 (E a /A ), (4) 
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where Ao = Tlvf/R is the single-particle level spacing 
of the dot and ugN labels different solutions for given 
(m,T, s). Equation (j4|) is easily solved numerically and 
the eigenstates to energy E a > are 



A a e 



imcp 



(5) 



where k a = E a /hvF and the normalization factor is 

A a = [x(Jm — Jm-lJm+1 + Jm+1 ~ JmJm+s)] (6) 

with J m = J m {E a / Ao). Time-reversal invariance 
implies the Kramers degeneracy relation E n m T S = 
E n _ m _i_ T _ s . Negative-energy (hole) solutions, 

<I>~ (r, <f>), follow by using the electron- hole symmetry 
property of the Hamiltonian, £7 n ,m,-T,s = —E n ^ m , T _ s . We 
use the multi-index a (a) to count states with positive 
(negative) energy. There is no zero-energy solution for 
B = 0, and we have a finite gap around the Dirac point . 

Next we add the magnetic field. Expressed in terms 
of the eigenstates < &i + ' ) and <3?~ \ the vector potential 
part in Hq has a matrix structure diagonal both in the 
quantum numbers (m, r, s) and the conduction/ valence 
band index ±, i.e., only different n states are mixed. By 
numerical diagonalization, it is straightforward to obtain 
the resulting eigenenergies E a > and Ea < 0, and the 
corresponding eigenstates. The indices n and thus a (a) 
are redefined to take into account the unitary transfor- 
mation diagonalizing Hq- Finally, we include the Zee- 
man term. Choosing the spin quantization axis along 
B, where s = ±1 corresponds to spin- up or spin-down 
states, the full eigenenergy E a > is given by^ii 



E„ — E n 



sfi B B, 



(7) 



and similary for E a < 0. In a slight abuse of notation, 
E a now denotes the full eigenenergy and not the solu- 
tion to Eq. ((4]) anymore. The Zeeman term is generally 
quite small^ but breaks the spin degeneracy of the levels, 
while the vector potential breaks the valley degeneracy, 
see Eq. (|3|). The resulting eigenstates are denoted by 
$i +) (r,0) and 



B. Many-body interactions 

We now include the Coulomb interaction among the 
particles. The nonintcracting reference problem is char- 
acterized by a filled Dirac sea (fi = 0), i.e., all E a < 
states are filled. The QED Hamiltonian H describing this 
problem can be expressed in terms of electron annhilation 
operators, c a , corresponding to the single-particle states 
3>o , and hole creation operators, <i~, with single-particle 
states <E>~ ^ . The full field operator is written as 

*(r) = E *i+'(r)c a + £ #|f >(r)4 (8) 

a a 

Since the Hamiltonian commutes with t z and s ■ B, we 
can write ^(r) = ^ rs ^ TS (r). The Hamiltonian is then 



given by H = Hk + Hi, with the kinetic part (note that 
E a <0) 



H k = J2E a clc a +J2\E a \4d a 



and the interaction term 

Hvpa ^ f drdr' 



2 ' ./ |r — r' 

tt' ss' 



(9) 



(10) 



where the colons denote normal ordering. Inserting the 
field operator expansion ([8]) into Eq. (fT0|) . 



H — i/fix + H 



(11) 



-fffix commutes separately with both the electron and the 
hole number operator, N e = J2 a c a c a an d Nh = Jia 4^a- 
The full Hamiltonian, however, only commutes with 
N = N e — Nh- We thus define the TV-particle problem by 
having A" excess electrons and N e h = Nh electron-hole 
pairs on top of the filled Dirac sea. Only N is conserved, 
while N e h can fluctuate. Under alone, the number 
iV e /j of electron-hole pairs is conserved, 

-fffix = H k + - E {Vabb'a' ~ 5ss'V aba ' ')c\c\c b ,C a , 
aba'b' 

+ \ E ^aWa'-^s'y aia r h )d\d\d y d a , (12) 
aba'b' 

- E ( V aWa' - S ss'V aba , h ,)cidld i ,C a ,. 
aa' ,bb' 

All terms not commuting with N e ^ are collected in the 
remaining part, H' = h + , with 

h = \ E (^ S '-*«'Km'&)4444 ( 13 ) 

aba'b' 

+ O E S ^-' V aba'b' C i d l d b C b' 



'bb' 



E {Vabb'a' - 5ss'V b ab'a')c a d\ 



twt J, 



- b <-b' 



abb' a' 



aa'bb' 



In Eqs. (|12j) and (|13[) . the spin quantum numbers are 
given by s = s a = s a > and s' = s b = s b > (when hole states 
are involved, a — > a etc.). These spin selection rules 
are encoded in the interaction matrix elements Vaa'b'b 
which have been derived in a form useful for numerical 
evaluation in Ref. [2^. We quote them for the convenience 
of the reader next. A finite matrix element follows only 
when the valley selection rule, r a = T a > and r b = Ty, and 
angular momentum conservation, m a + m a i = m b + m b > , 
are satisfied. When all selection rules are met, 
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aa'b'b 



= (4ir) 2 aA A a A a ,A b ,A b ^2c q ,i / dr r - ' (J ma (E a r) J mb {E b r) + J ma+x {E a r) J mb+1 (E b r)) 

1=0 Jo 

x f dr' (r') l+1 (J ma , (E a ,r')J mb , {E v r>) + J ma , +1 (E a ,r')J mb , +1 (E b ,r')) , 
Jo 



with E a in units of Ao = Hvf/R and q = \m b — m a \. The 
coefficient Cq i vanishes when I + q is odd or when I < q. 
For q = I = 0, we have Co o = 1/2, while otherwise 

_ (2/-l)!! ( '^ /2 (n-l/2)(n-l-l) 
~ 2 l +H\ 11 



n{n~l- 1/2) 



C. Calculation approaches 

A standard way to proceed is to employ the no-pair 
approach^ With the projector A + to the subsector E a > 
of the single-particle Hilbert space (for each particle), 
we thus consider the TV-particle problem with respect to 
the filled Dirac sea. The projected Hamiltonian H + = 
A + HA + = A + Hfi x A + is given by 

H+ = £25,4^ (14) 

a 

+ ^ ^2 ( V abb'a' ~ S ss <Vaba>b')clclc b ,C a ,. 
aba'b' 

As detailed in Refs. I25TI26L this allows for a straightfor- 
ward implementation of the HF approach, and we will 
report HF results in Sec. IIVI In contrast to Refs. [25ll26l . 
we here include the valley and spin degrees of freedom. 
Given the converged self-consistent density matrix, one 
can obtain the ground-state energy E(N), the total spin 
quantum number S of the TV-electron dot from the eigen- 
values h 2 S(S+ 1) of YliLi S?, and the valley polarization 
eigenvalue t(N) = J2i T i- 

On a more general level, we allow for a fixed number 
of electron-hole pairs by employing the Hamiltonian Hj\ x 
[Eq. ([T2])]. The no-pair Hamiltonian H+ follows from 
with N e h = 0. When the ground-state energy of Hs x is 
minimized for some N ebi > 0, the no-pair approach breaks 
down. Interactions are then able to overcome the gap 
between valence and conduction band, and one cannot 
treat the Dirac sea as inert anymore. Hfi x as well as 
the QED model H are studied by ED and perturbation 
theory in Sec. IIIII 



III. PARTICLE-HOLE PAIR PRODUCTION 
AND RECONSTRUCTION OF THE GROUND 
STATE 

We now compare the three theoretical approaches in 
Sec. Ill CI by employing ED for particle numbers N = 2 




FIG. 1: (Color online) ED results for the ground-state energy 
E in units of Ao = Fivf/R vs fine structure constant a for 
N = 2 particles in a graphene dot with 5 = 0. We here use 
the Hamiltonian fig* in Eq. (|12[) , 



and 3. Convergence in the ED calculations was achieved 
by keeping about 30 single-particle states (per spin and 
valley degree of freedom), and memory size limitations 
represented the main bottleneck. For a given a, ED re- 
sults can be obtained within a few minutes on a standard 
desktop computer. 

Figure [T] shows results for the a-dependence of the 
ground-state energy for N = 2 using H& x with N e h = 
and 1. We observe that for a < 1.3, the ground state 
of -f/fix contains no electron-hole pair, but for stronger 
interaction, the ground state undergoes reconstruction 
and involves at least one electron-hole pair. The no-pair 
Hamiltonian H+ thus necessarily fails when a > 1.3. 
Moreover, as we shall discuss next, the presence of H' 
[Eq. ([i"3]) ] restricts its applicability even further. 

Since ED of the QED Hamiltonian H in Eq. (p} is 
computationally very demanding even for N = 2, in the 
remainder of this section, we shall restrict ourselves to a 
spinlcss single- valley version of graphene. The ED results 
obtained from 77fi x and H are compared for N = 2 in 
Fig. [21 For the spinless single- valley version of Hn x , no 
electron- hole pairs are excited in the ground state for a < 
1.9. However, the H' contribution is important already 
for a > 0.5, see Fig. [2j The full interaction correction to 
the energy is significantly lowered by including H' and 
may even change sign for large a. In these calculations, 
the Hilbert space was truncated to contain at most one 
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FIG. 2: (Color online) ED results for the interaction energy 
5E(a) = E(a) - E(0) (in units of A ) vs a for N = 2. 
We consider a spinless single-valley version of graphene 
with B = 0. The curves for iV e h = 0, 1, 2 correspond to 
the Hamiltonian _fffi x with N e h electron-hole pairs, i.e., the 
ground state then has no electron-hole pair for a < 1.9. 
However, the full QED Hamiltonian (|11[) . where we truncate 
the Hilbert space to at most one or two electron-hole pairs 
(max(iV e h) = 1,2), has a significantly lower energy already 
for a > 0.5. 
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FIG. 3: (Color online) Same as Fig. [2] but for N = 3. 



or two electron-hole pairs. For a < 1.5, this appears to 
be sufficient. Figure [3] shows results for N ~ 3, where we 
arrive at similar conclusions. 

The effect of H' can also be evaluated analytically by 
using second-order perturbation theory (the first order 
vanishes identically). The result is shown for N = 2 in 
Fig. |4l together with the ED results from Fig. [2j We 
see that second-order perturbation theory captures the 
ED data quite well, especially for a < 1. The same 
conclusion was reached for N = 3 (results not shown 
here), and the combination of ED (or HF) calculations 
for _fffi x supplemented with a perturbative treatment of 
H' should in general provide a good approximation of 
the ground state. 
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FIG. 4: (Color online) Same as Fig. [2] but including the 
results of second-order perturbation theory in H' . The 
curve AE\ takes into account corrections involving one 
electron-hole pair only, while AE12 is the full second-order 
result including up to two pairs. 
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FIG. 5: (Color online) Main panel: Relative number of 
electron-hole pairs in the ground state, {N e h)/N, vs a for 
N — 2 and several values of the magnetic field. We take the 
dot radius R = 30 nm. The results were obtained using the 
QED Hamiltonian but with the Hilbert space truncated to 
at most two electron-hole pairs. Inset: Interaction correction 
SE vs a for N = 2 and B = 1 T. Shown are ED results using 
H+ (N e h = 0) and for the full H, where the Hilbert space 
was truncated at max(Af e t) = 2. 



Let us now discuss the case of finite magnetic field, 
again for the computationally simpler spinless single- 
valley case with N = 2. We have also studied N = 3 
particles, again with very similar results. The main panel 
of Fig. [5] shows the average number of electron-hole pair 
excitations in the ground state for several values of the 
magnetic field. The shown results are for a dot radius 
R = 30 nm. A magnetic field of B = 1 T corresponds 
to the magnetic length Ib = (c/eB) 1 / 2 ss 26 nm, which 
is of the same order of magnitude as the radius. Evi- 
dently, in a magnetic field, the proliferation of electron- 
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FIG. 6: (Color online) HF results using the no-pair model, 
H+, for the ground-state spin S vs particle number N for 
B = and for B = B z = 3 T (with R = 30 nm). For all 
shown JV and a < 1, S(N) does not depend on a. 




FIG. 7: (Color online) HF results for the ground-state valley 
polarization t(N) = J^fLi Ti m the zero-field case with a — 
and a = 1. For a = 0.5, the same result as for a = 1 is 
found. 



hole pairs becomes more important. We interpret this 
effect as the finite-size analogue of the magnetic catalysis 
phenomenon^ The interaction correction to the ground- 
state energy is shown for B = 1 T in the inset of Fig. [5] 
While the result shows qualitatively similar behavior as 
for B = 0, the now more significant deviations between 
the ED data and the no-pair result are consistent with 
magnetic catalysis again. We note in passing that for 
a > 1.5, the basis size used in our ED calculations is 
most likely not sufficient, and probably N e h > 2 states 
also contribute to the ground state. The stcplike features 
in Fig. [5] are then presumably smeared out. 

We conclude that the no-pair Hamiltonian is quantita- 
tively reliable only for weak interactions, a < 0.5, and for 
not too large magnetic fields. For stronger interactions 
and/or fields, the ground state undergoes reconstruction 
and electron-hole pair proliferation. Using [Eq. (fll2"j)] 
is not sufficient to get more accurate results, and one has 
to include H' [Eq. (fj~3|) ] which does not conserve the elec- 
tron and hole numbers separately. However, for a < 1, 
quite accurate results for the ground-state energy are ob- 
tained by combining ED (or HF) calculations for the no- 
pair Hamiltonian with subsequent second-order pertur- 
bation theory in H' . Here only two or three particles 
have been addressed, where electron-hole proliferation 
takes place around a k. 1. Since the bulk case, which 
follows from the above model by a suitable limiting pro- 
cedure with N — > oo and R — > oo, has a phase transition 
at a ~ 1.1, the finite-size crossover apparently depends 
on N only weakly. 



IV. ADDITION SPECTRUM AND 
GROUND-STATE PROPERTIES 

Let us now turn to HF results for the ground state 
of the A^-electron dot with N < 15 using the no-pair 
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FIG. 8: (Color online) Same as Fig.[7|but for a perpendicular 
magnetic field with B = 3 T (with R = 30 nm). 
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FIG. 9: (Color online) Addition energy (|15[) in units of Ao 
vs N for several a from HF calculations for H+ with B — 0. 
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FIG. 10: (Color online) Same as Fig. M but with B x = 15 T 
and B Z = 3T (for i? = 30 nm). 



Hamiltonian i7 + [Eq. (|H|) ]. As discussed in Sec. IIII1 this 
approximation is reliable only for weak interactions, and 
we focus on the regime a < 1 below. The spin and valley 
degrees of freedom are fully included in our self-consistent 
HF calculations. 

The total ground-state spin S follows from the eigen- 
value h 2 S(S + 1) of the total squared spin operator and 
is shown as a function of N in Fig. [51 both for B = 
and in the perpendicular magnetic field B = 3 T. For 
N < 14 and a < 1, the spin filling sequence S(N) is 
independent of the interaction strength a and displays a 
four-periodicity for B = 0. For B ^ 0, this periodicity 
is reduced to a two-periodicity since now spin degener- 
acy is broken. We note that the spin filling sequence 
can be measured experimentally by Coulomb blockade 
spectroscopy^ 

The total ground-state valley polarization, t(N), is 
shown in Fig.[7]for B = 0, and in Fig.[8]under the perpen- 
dicular field B = 3T. When B = 0, the full Hamiltonian 
is symmetric under r — > — r, and we here show only the 
positive solution. However, a finite orbital field B lifts 
the valley degeneracy. We observe from Fig. [7] that for 
B = 0, interactions reduce the four-periodicity of t(N) 
for a = down to a two-periodicity. This can be un- 
derstood by noting that the interaction of particles in 
different valleys is typically weaker than the intra- valley 
interaction. For B ^ 0, this implies pronounced inter- 
action effects on the valley polarization. Figure O shows 
that strong interactions force subsequent particles to be 
added into the same valley, thereby valley-polarizing the 
iV-particle system. Both the spin and valley filling se- 
quences obtained by HF theory have been independently 
confirmed by ED of the no-pair Hamiltonian for N < 4 
(data not shown). 

To estimate the accuracy of the HF approximation for 
the no-pair model, we have also determined the relative 
difference between the HF (£hf) and the ED (-Eed) cn- 



Ej3f(N, a)-E ED (N,a) 
E EB (N 1 a)-E EU (N,0) 



In all studied cases (N < 4), S(N) was found to be rather 
small. For instance, even when taking the large value 
a = 1.5, we obtain 8(2) = 0.107, 5(3) = 0.175 and 6(4) = 
0.148. As long as the no-pair approach stays valid, we 
conclude that HF theory yields quite accurate results. 

Our HF results for the addition energy^ which follows 
from the ground-state energy E(N) using the relation 

A(iV) = E(N + 1) + E(N — 1) — 2E(N), (15) 

are shown in Fig. [9] for B = and several a. (Simi- 
lar HF results but for the spinless single-channel version 
were discussed in Ref. H^.) Peaks in A (AT) signify espe- 
cially stable dot configurations (magic numbers). While 
for a = 0, A(A) again shows the four-periodicity due 
to spin- valley degeneracy, the addition energy peaks be- 
come less pronounced with interactions, and the four- 
periodicity is not always visible. Interestingly, while 
there are magic numbers N = 4, 8, 12, . . . related to com- 
pletely filled "energy shells" in the noninteracting case, 
the addition energy curves A (A) are rather featureless 
and almost flat for strong interactions. This indicates 
that a constant interaction models provides a reasonable 
description, where the microscopic Coulomb interaction 
is replaced by the electrostatic charging energy of an ef- 
fective capacitor. The addition energy A(A) for B ^ is 
shown in Fig. 1101 The in-plane part B x of the magnetic 
field here acts to increase the spin Zeeman field. How- 
ever, Zeeman effects in graphene are weak,— and indeed 
almost the same results as those in Fig. [TU] were found for 
B x — and B z = 3 T. As a consequence of the broken 
spin degeneracy, only an (approximate) two-periodicity 
in A(A) is observed in the magnetic field case. 



V. CONCLUSIONS 

In this work, we have studied the ground state prop- 
erties of N particles in a disorder-free circular graphene 
quantum dot, with the filled Dirac sea as the point of 
reference. The boundary of the dot has been mod- 
elled by the infinite-mass boundary condition, and the 
particles interact via the unscreened Coulomb potential 
whose prefactor is proportional to the bare dimensionless 
fine structure constant a. In contrast to atomic physics 
where a = 1/137 is very small, in graphene (e.g., by the 
variation of the substrate dielectric parameter) a may 
be tuned up to a maximum value of a ~ 2.2 (reached 
for freely suspended samples). For instance, a recent 
experiment^ using Coulomb blockade spectroscopy for 
a graphene dot reported awl. We have studied the 
A-particle problem in a graphene dot on various levels 
of complexity - from the no-pair Hamiltonian to the full 
QED model - and by a number of different techniques. 
Our main results arc as follows. 
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By using exact diagonalization (ED) for TV = 2 and 
3 particles, we found that the no-pair Hamiltonian H + 
originally proposed by Sucher,— where the filled Dirac 
sea is assumed to be inert, is quantitatively reliable only 
for a <C 1 , see Sec. IIIIl While this represents the stan- 
dard situation in atomic physics,— it can easily be vi- 
olated in graphene. For a > 0.5, our calculations in- 
dicate that electron-hole pair excitations contribute to 
the ground state energy. For a > 1, these excitations 
proliferate and eventually cause a completely restruc- 
tured ground state. Technically, the projection opera- 
tor A + defining the vacuum should thus be changed to 
include interaction effects in a self-consistent manner. 
Mittleman^ has shown that this goal can be achieved 
by first minimizing the ground state energy E(N,A+) 
for given A + , followed by the maximization of the energy 
over all possible A+. The final result for E(N) should 
then be equivalent to the QED results obtained numeri- 
cally by ED (in the limit of infinite basis size) . 

We here argue that graphene dots realize a finite-size 
crossover version of the bulk semimctal-insulator phase 
transition. We find that the crossover scale is set by 
awl, consistent with the bulk result a c « 1.1 £ When an 
orbital magnetic field is applied - the Zeeman field plays 
no significant role - electron-hole pair proliferation sets 
in earlier and implies a lowering of a c , consistent with 
the magnetic catalysis scenario^ Even on a qualitative 
level, the no-pair Hamiltonian H + is thus reliable only 
on the semimetallic side of the transition (a < 1). 

For the regime a < 1, we have reported detailed results 
using Hartree-Fock theory for H + and N < 15 particles 
in Sec. IIV1 taking into account the spin and valley de- 



grees of freedom. We find a four- (two-)periodicity in 
the spin filling sequence in the absence (presence) of a 
magnetic field, which can be understood from the single- 
particle picture and remains unaffected by weak inter- 
actions. However, the valley filling sequence is more in- 
tricate, especially when B / 0. This is related to sub- 
tle differences between the intra- and inter-valley scat- 
tering matrix elements of the Coulomb interaction. We 
observe a strong tendency towards valley polarization in- 
duced by interactions in this A^-body problem. Finally, 
our analysis of the addition energy spectrum reveals that 
the constant interaction model can provide a reasonable 
description. 

In our previous HF study of the spinless single-valley 
no-pair problem^ we found that Wigncr molecule for- 
mation sets in for strong interactions. Since that regime 
corresponds precisely to the onset of electron-hole prolif- 
eration, a > 1, where the no- pair model becomes unreli- 
able, we have analyzed the question of Wigner molecule 
formation using ED for A~ = 3 under the full QED 
model again. The Wigner molecule is identified from pro- 
nounced density correlations, and our numerical results 
(not shown here) are very similar to what we reported in 
Ref. We thus expect that the Wigner molecule for- 
mation is only weakly affected by the electron-hole pair 
proliferation reported in this paper. 
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